Image registration method

ABSTRACT

Disclosed is an image registration method by iteratively determining a transformation that is optimal regarding a given distance criterion and smoothness criterion. Said method allows corresponding landmarks in the images to be definitely represented on top of each other and comprises the following steps: (1) an iteration counter and the initial displacement field are initialized; (2) the numerical solutions of the non-linear partial differential equation (PDE) are determined by means of the differential operator that can be derived from a predefined smoothness criterion and the point evaluation functionals located at given landmarks; (3) the interpoliation conditions are combined; (4) a special numerical solution of the PDE is calculated by means of the force that is determined based on the distance criterion and the actual displacement field as well as the differential operator derived from the smoothness criterion; (5) the special solution is evaluated at the landmarks; (6) the coefficients are calculated in order to calculate an updated displacement; (7) the displacement field is updated and the iteration counter is increased; (8) the displacement is verified regarding convergence; and (9) steps (4) to (8) are repeated if the convergence criterion is not satisfied.

The invention relates to an image registration or recording method, i.e.for the correction of geometrical differences in differentrepresentations of an object. These methods play an important part, e.g.in medical technology and particularly when analyzing tissue changes inconjunction with early cancer diagnosis.

Methods are already known which carry out an image registration on thebasis of a distance criterion (Lisa Gottesfeld Brown: A survey of imageregistration techniques, ACM Computing Surveys, 24(4): 325-376, 1992,Jan Modersitzki: Numerical Methods for Image Registration, Habilitation,Institute of Mathematics, University of Lubeck, Germany, 2002). Thegeneral methodology is based on the optimization of a target function tobe chosen in use-conforming manner and which is typically based on imageintensities. In such methods, apart from the image information, nofurther information is used for registration purposes. The registrationresult is only of an optimum nature in the sense of a global averaging.If particular significance is attached to specific, characteristicpoints in an application (such as e.g. the so-called anatomicallandmarks in medical applications) such methods cannot be recommended.

Apart from image registration on the basis of a distance criterion,methods are also known which perform image registration exclusively onthe basis of control points (Karl Rohr: Landmark-based Image Analysis.Computational Imaging and Vision. Kluwer Academic Publishers, Dordrecht,2001). In such methods prospectively or retrospectively correspondingcontrol points are associated with the views to be registered and arethen matched by means of registration. The disadvantage of such methodsis that the registration exclusively takes account of control points.Such methods cannot take account of further image informations, such ase.g. image intensities. In the case of unsatisfactory registrationresults a user can only attempt to improve them by skilled introductionof further control points. The insertion of further control points isbased on subjective trial and error for which no guidelines exist and inparticular there is no automated procedure.

The problem of the invention is to develop an image registration method,which leads both to a perfect, guaranteeable correspondence between anumber of predetermined control points and also an optimum result in thesense of the distance criterion.

According to the invention this problem is solved by the iterativedetermination of a transformation optimum with respect to apredetermined distance and smoothness criterion, in which control pointscorresponding in the images are imaged on one another in guaranteeablemanner by (1) initializing an iteration counter and the initialdisplacement, (2) determining the numerical solutions of the nonlinear,partial differential equation (PDE) with the differential operatorderivable from a predetermined smoothness criterion and the pointevaluation functionals located at the predetermined control points, (3)combining the interpolation conditions, (4) calculating a specific,numeral solution of the PDE with the force determined on the basis ofthe distance criterion and the actual displacement field and thedifferential operator derived from the smoothness criterion, (5)evaluating the specific solution at the control points, (6) determiningthe coefficient for calculating an updated displacement, (7) updatingthe displacement field and raising the iteration counter, (8) checkingthe displacement for convergence and (9) in the case of nonfulfilment ofthe convergence criterion repetition of steps (4) to (8).

The method sequence is illustrated by the flow chart of FIG. 1.

For simplification purposes a view is referred to as a reference image(reference R) and a further view, which is to be corrected, as atemplate (template T). From the formal standpoint these are functions ofa d-dimensional, real space or a subset Ω⊂R^(d) in the set of realnumbers. Thus, to each d-dimensional point xεQ is associated throughR(x) and T(x) a value which can be interpreted e.g. as a colour or greyvalue.

In practical applications, particularly during every programming of thepresent method, the reference and template can be in discreet form. Theimages are then functions on a lattice (e.g. Ω={1, . . . , n1}×{1, . . ., n2} for the dimension d=2) in a discreet set (e.g. in the set {0, . .. , 255}) and can be interpreted as being formed from pixels. For theregistration method these restrictions and in particular the specificnature of the discretization are unimportant. The restrictions aresolely made for simplified description purposes. The method can be usedin the same way on random d-dimensional data sets.

The function of image registration consists of the determination of adisplacement function u, so that the requirement R(x)=Tu(x) is optimumwell fulfilled with the short form Tu(x):=T(x−u(x)) for all xεΩ. Forcalculating the template Tu deformed by u in the case of discreet,predetermined images such as are of a conventional nature in imageprocessing an interpolation (e.g. d-linear) has to be performed, becausethe displaced coordinates x−u(x) are not necessarily located on thediscreet lattice. The way in which such an interpolation takes place isunimportant for the registration method.

Over and beyond the aforementioned similarity, requirements must be madeon the displacement smoothness and on the imaging characteristics withrespect to a number of preselected control points. In the simplest casethe coordinates of each of the m control points K^(T.j) of the templatemust be imaged on the in each case corresponding control point k^(R.j)of the reference, j=1, . . . , m. If the coordinates of the controlpoints coincide, which may be ensurable by a preregistration, then u=0applies at these points.

As is normally the case with optimization problems, the determination ofa minimizer of the aforementioned distance criterion can take placeiteratively by means of a gradient descent method. In principle, anyrandom distance criterion can be selected. The forces associated withthe standard distance criteria are given in the literature (Modersitzki2002). The specific nature of the calculation of these forces isunimportant for the registration method.

In principle, any functional known from the literature can be used asthe smoothness criterion. From the smoothness criterion it is possibleto derive a partial differential operator A. These operators are knownfor the criteria used in the literature (Modersitzki 2002). The soughtdisplacement u can then be characterized as a solution of a nonlinear,partial differential equation (PDE).

For determining a numerical solution of this PDE use is made of a finitedifferential approximation of the differential operator, which thenleads to an equation system for the lattice values of the displacement.However, the specific discretization of the differential equation lackssignificance for the registration method.

This procedure coincides with the method based solely on the distancecriterion and the smoothness criterion. The new aspect consists of asuitable binding in of the predetermined control points into thedisplacement calculation, in which a correspondence of the controlpoints can be guaranteed. As methods are already known for thedetermination of the displacement based on the distance and smoothnesscriterion, a method is given here which combines partial solutions in anappropriate manner so as to give an overall solution, e.g. in the form${{{ul}(x)} = {{v^{o}{l(x)}} + {\sum\limits_{j = 1}^{m}{\lambda_{j}^{l}v^{j}{l.(x)}}}}},\quad{x \in \Omega},\quad{l = 1},\ldots\quad,d$

A is the differential operator associated with the smoothness term and fis the field of forces belonging to the distance criterion, so thatv^(o) is a numerical solution of Av^(o)=−f, the functions v^(j) arenumerical solutions of the distributional PDE Av^(j)=δ^(j), j=1, . . . ,m, in which δj locates the point selection functional (Dirac impulse) atcontrol point K^(T.j). The specific, numerical method for the solutionof the PDE is unimportant for the registration method.

From the mathematical standpoint v^(j), j=1, . . . , m are Green'sfunctions of the differential operator A, representing a solution of thePDE at the predetermined single point displacement. A suitable linearcombination of these Green's functions consequently ensures that in theoverall method, in the required manner, all the control points areimaged on one another.

The function v^(o) is so determined by means of an iterative method thatthe distance criterion is minimized, whilst maintaining the requiredsmoothness. The weight functions λ^(l)j are so adapted that the controlpoints are imaged in the required way.

The initialization of the program requires the selection of a distanceand a smoothness criterion or the force derivable from said criteria andthe differential operator. On the basis of the point evaluationfunctionals located at the control points it is then possible todetermine the Green's functions v^(j), j=1. m with a numerical method.They are not changed during the further course of the method.

The inventive initialization is followed by a standard iterationprocedure, during which there is a gradient descent, whilst takingaccount of the control points. No human intervention is needed. Thus,the described method combines the advantages of methods based ondistance criteria (particularly automatability and an on average optimumregistration) with those of the control point method (guaranteedregistration of distinguished points) and on predetermining an initialset of control points gives reproducible, optimum results, independentof the user or computer program. No important part is played by computercode details in the final result of the image registration and they onlyinfluence the requisite computing time and memory requirements.

The images to be registered can be digital images, pixels, JPEG,wavelet-based objects or acoustic signals.

The linear equation systems occurring in the method can be solveddirectly, indirectly, iteratively or by multigrid and for the method usecan be made of a reference coordinate system imaged by Euler or Lagrangecoordinates.

The invention also proposes the registration of one, two orthree-dimensional and sequences of one, two and three-dimensionalobjects, as well as the use of control points in the form of anatomicallandmarks, fiducial markers or other characteristic quantities.

The distance criterion proposed is based on intensity, edge, corner,surface normal or level set or on the “sum of squared differences”, L2distance, correlation, correlation variants, mutual information ormutual information variants.

The force terms associated with the distance quantity are to becalculated by finite difference methods or gradient formation and thesmoothness criterion used is to be physically motivated by means of anelastic potential or a fluid approach or diffusive or curvatureapproaches based on time or space derivatives of the displacement.

The boundary conditions of the differential operator are advantageouslygiven by explicit or implicit, Neumann, Dirichlet, sliding, bending orperiodic boundary conditions.

The nature of the discretization of the differential operator should bebased on finite differences, finite volume, finite elements, Fouriermethods, series expansions, filter techniques, collocations or multigridand interpolation is to be performed d-dimensionally by means of splinesor wavelets.

Finally, displacement can be explicitly updated by means of theincrement of the displacement or its time derivative.

1. Method for the registration of images by iterative determination ofan optimum transformation with respect to a predetermined distance andsmoothness criterion, characterized in that control points correspondingin the images can be imaged on one another in guaranteeable manner, by(1) initialization of an iteration counter and the initial displacementfield, (2) determining the numeral solutions of the nonlinear, partialdifferential equation (PDE) with the differential operator derivablefrom a predetermined smoothness criterion and the point evaluationfunctionals located at the predetermined control points, (3) combiningthe interpolation conditions, (4) calculating a specific, numericalsolution of the PDE with the force determined on the basis of thedistance criterion and the actual displacement field and thedifferential operator derived from the smoothness criterion, (5)evaluating the specific solution at the control points, (6) determiningthe coefficients for calculating an updated displacement, (7) updatingthe displacement field and raising the iteration counter, (8) checkingthe displacement for convergence and (9) in the case of nonfulfilment ofthe convergence criterion, repetition of steps (4) to (8).
 2. Methodaccording to claim 1, characterized in that one, two orthree-dimensional objects, as well as sequences of one, two andthree-dimensional objects are registered.
 3. Method according to one ofthe preceding claims, characterized in that the control points areanatomical landmarks, fiducial markers or other characteristicquantities.
 4. Method according to one of the preceding claims,characterized in that the distance criterion is based on intensity,edge, corner, surface normal or level set or on the sum of squaredifferences, L2 distance, correlation, correlation variants, mutualinformation or mutual information variants.
 5. Method according to oneof the preceding claims, characterized in that the force termsassociated with the distance quantity are calculated by means of finitedifference methods or gradient formation.
 6. Method according to one ofthe preceding claims, characterized in that the smoothness criterionused is physically motivated by means of an elastic potential or a fluidapproach or on diffusive or curvature approaches based on time or spacederivatives.
 7. Method according to one of the preceding claims,characterized in that the boundary conditions of the differentialoperator are explicit or implicit, Neumann, Dirichlet, sliding, bendingor periodic boundary conditions.
 8. Method according to one of thepreceding claims, characterized in that the nature of the discretizationof the differential operator is based on finite differences, finitevolume, finite elements, Fourier methods, series expansions, filtertechniques, collocations or multigrid.
 9. Method according to one of thepreceding claims, characterized in that the interpolation is performedd-dimensionally by means of splines or wavelets.
 10. Method according toone of the preceding claims, characterized in that the displacement isexplicitly updated by means of the increment of the displacement or itstime derivative.